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NUMERICAL  SIMULATION  OF  THE  COMPRESSIBLE  ORSZAG-TANG 
VORTEX  D.  SUPERSONIC  FLOW 

I.  INTRODUCTION 

In  a  recent  paper,  to  be  designated  as  I  [1],  we  initiated  a  systematic  research  program 
Ubiiig  numerical  simulations  of  dissipative,  fully  compressible  magnetofluids  to  identify  and 
analyze  the  effects  of  compressibility  on  magnetohydrodynamic  (MHD)  turbulence.  The 
first  studies,  reported  therein,  tracked  the  decay  of  initially  subsonic  flow  fields  and  used 
the  analogous  incompressible  flows  as  baselines  for  comparison.  The  initial  average  Mach 
number  M  ranged  from  0.2  to  0.6,  and  the  ratio  of  thermal  to  magnetic  pressure  (/3)  ranged 
from  30.0  to  10/3.  The  initial,  dimensionless  mass  density  was  uniform  at  p  =  1  while  the 
Lundquist  numbers  varied  from  50  to  200. 

The  present  paper  deals  with  the  effect  of  embedded  supersonic  flows  and  the  re¬ 
sulting  emerging  shock  waves  on  phenomena  associated  with  MHD  turbulence,  including 
reconnection,  the  formation  of  current  sheets  and  vortex  structures,  and  the  evolution  of 
spatial  and  temporal  correlations  among  physical  variables.  The  parameters  used  here  are 
M  =  1.0  and  1.5,  p  =  2.78  and  6.25,  respectively,  and  =  10/3,  with  the  same  range  of 
Lundquist  numbers  as  above.  Often  flows  within  regions  bounded  by  shocks  are  locally 
subsonic,  and  the  previous  subsonic  calculations  should  apply  to  those  regions.  How¬ 
ever,  both  the  supersonic  disturbances  responsible  for  forming  the  shocks  and  the  current 
sheets  accompanying  the  shocks  will  have  effects  on  the  overall  flow  field  which  cannot 
be  predicted  on  the  basis  of  the  subsonic  studies  alone.  An  example  is  the  production 
of  turbulence  by  inhomogeneous  laser  pulses  propagating  in  a  gas  [2,  3].  The  laser  pulse 
produces  a  channel  in  the  background  gas,  with  some  regions  of  the  channel  cross  section 
heated  more  than  others.  The  resulting  shocks  interact  with  the  inhomogeneous  density 
distribution  of  the  channel  to  produce  turbulence.  This  paper  provides  initial  descriptions 
and  assessments  of  similar  supersonic  effects  upon  a  compressible  magnetofluid.  Our  re¬ 
sults  should  be  useful  in  describing  local  supersonic  regions  in  a  turbulent,  compressible, 
dissipative  magnetofluid. 

Among  the  obvious  complicating  factors  in  fully  compressible  MHD  flows  are  density 
and  entropy  fluctuations,  finite  magnetoacoustic  waves,  discontinuities  (e.g.,  shocks  and 
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contact  surfaces),  and  a  finite  information  transport  velocity.  Theories  of  “weakly  com¬ 
pressible”  MHD  flows  [4,  5,  6,  7],  while  both  useful  and  elegant,  account  only  for  density 
fluctuations  (and  not  entropy  fluctuations),  filtering  acoustic  waves  out  of  the  solution  or 
damping  them  out  on  a  rapid  time  scale.  To  reveal  the  full  range  of  compressible,  turbulent 
phenomena,  initial  studies  must  necessarily  survey  the  phenomenology  of  relevant  nonisen- 
tropic  flows  which  include  all  of  the  effects  mentioned  at  the  beginning  of  this  paragraph. 
One  can  then  begin  to  identify  problems  requiring  theoretical  treatment  and  exhaustive 
numerical  study  using  the  maximum  achievable  resolution. 

Our  approach  differs  from  that  of  many  other  studies  of  turbulence  in  that  we  have 
chosen  a  model  problem,  the  Orszag-Tang  vortex  system  [8],  which  involves  decay  from 
nonrandom  initial  conditions,  rather  than  to  follow  the  decay  of  a  random  initial  field  or 
the  evolution  of  a  driven  turbulent  flow.  The  Orszag-Tang  system  is  doubly  periodic,  and 
the  initial  conditions  consist  of  single-mode,  solenoidal  velocity  and  magnetic  fields,  each 
containing  X-points  and  O-points.  The  initial  mass  density  is  flat  and  the  initial  pressure 
fluctuations  are  incompressible,  balancing  the  local  forces  for  a  magnetofluid  of  unit  mass 
density.  The  identification  and  tracking  of  the  various  compressible  MHD  phenomena 
thus  are  more  straightforward  than  in  a  statistical  approach.  In  addition,  such  a  simple, 
easily  reproducible,  model  problem  can  serve  as  a  useful  reference  for  comparison  with 
the  results  from  studies  of  more  realistic  turbulent  fields  [9].  This  model  problem,  in 
fact,  contains  most  of  the  significant  features  of  MHD  turbulence,  including  dissipation  of 
magnetic  and  kinetic  energy,  reconnection,  formation  of  high-density  jets,  selective  decay, 
dynamic  alignment,  and  the  emergence  and  manifestations  of  small-scale  structure.  The 
incompressible  version  has  received  thorough  study  [10]  and  provides  an  excellent  baseline 
for  direct  quantitative  and  structural  comparisons  with  the  evolution  in  a  compressible 
medium.  Past  treatments  also  have  emphasized  the  relationship  of  the  Orszag-Tang  vortex 
system  to  the  theory  of  resistive  tearing  modes  [11].  Similar  highly  resolved  compressible 
calculations  should  be  performed  for  various  magnetoacoustic  Mach  numbers  and  values 
of  plasma  /?.  The  popularity  and  simplicity  of  the  model  also  make  it  an  obvious  choice 
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for  comparison  and  calibration  of  numerical  algorithms  for  simulating  compressible  MHD 
turbulence. 

Our  previous  paper  [1]  has  added  subsonic  flow  as  another  useful  baseline  for  studies 
of  supersonic  flows.  That  paper  provides  a  background  on  the  relevant  literature  and  a 
detailed  discussion  of  various  formulations  and  numerical  methods  used  for  compressible 
MHD  turbulence  and  associated  phenomena.  The  present  paper  deeds  with  situations  in 
which  regions  of  appreciable  supersonic  flow  lead  to  the  emergence  of  shocks.  These  flows 
have  a  major  effect  on  reconnection  at  the  X-points  in  the  magnetic  field,  the  density  jets 
originally  observed  at  subsonic  Mach  numbers,  the  formation  of  other  density  extrema, 
the  dissipation  of  energy  through  the  formation  of  current  sheets  at  the  shocks,  and  the 
spectral  composition  of  autocorrelations  and  cross-correlations  of  physical  variables.  The 
next  section  describes  briefly  our  formulation  and  numerical  methods.  Section  III  presents 
our  results  on  the  evolution  of  the  local  structure  of  the  flow  field,  the  global  properties 
of  the  system,  and  spectral  correlations.  Section  IV  discusses  the  important  dynamical 
properties  and  observational  consequences  of  embedded  supersonic  regions  and  emerging 
shocks  in  the  Orszag-Tang  model  of  an  MHD  system  undergoing  reconnection.  We  also 
draw  conclusions  regarding  the  effects  of  local  supersonic  regions  on  MHD  turbulence. 
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II.  NUMERICAL  SIMULATIONS 
A.  Formulation 

The  following  nonlinear  partial  differential  equations,  written  in  dimensionless  form, 
govern  the  behavior  of  a  two-dimensional,  compressible,  dissipative  magnetofluid  in  our 
numerical  simulations: 


dp 

dt 
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dt 

dB 

dt 
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dt 
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V  x 
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V-B  =  0, 


with  an  equation  of  state, 


p  =  (7-1)17  . 


(le) 


(2) 


Here  p  is  the  mass  density,  v  is  the  flow  velocity,  p  is  the  mechanical  pressure,  I  is  the  unit 
dyad,  B  is  the  magnetic  induction  field,  U  is  the  internal  energy  density,  E  is  the  total 
energy  density,  given  in  dimensionless  units  by 

E(x,t)m  p|v|2  +  |B|2+  V,  (3) 
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T  =  p/p  is  the  dimensionless  temperature,  r  =  [(djt;,  +  d{Vj)  —  |V  •  v£tJ]  e,-  ej  is  the 
stress  tensor,  and  7  =  5/3  is  the  ratio  of  specific  heats.  The  thermal  conductivity  (/c), 
magnetic  resistivity  (77),  and  viscosity  (p)  are  constant  and  uniform,  and  we  have  assumed 
that  the  bulk  viscosity  is  zero  [12].  Important  dimensionless  numbers  are  the  viscous 
Lundquist  number  Sv  =  poV^Lo/p  [13],  the  resistive  Lundquist  number  Sr  =  V%Lo/r), 
and  the  Prandtl  number  Pr  =  cpp/k. 

In  these  definitions,  po  is  a  characteristic  density,  V%  is  a  characteristic  Alfven  speed, 
Lo  is  a  characteristic  length  equal  to  the  reciprocal  of  the  minimum  wavenumber  kmm  = 
2n/L  (L  =  system  length),  cp  is  the  specific  heat  at  constant  pressure,  and  Vo  =  is  a 
characteristic  flow  speed.  The  thermodynamic  normalization  sets 

Eo=po=B20/&7r  =  poV¥/2  .  (4) 

Time  (t)  is  measured  in  units  of  the  Alfven  transit  time  Lo/V%.  For  the  runs  reported 
here,  we  set  Sv  =  Sr  =  50,  100,  or  200  and  Pr  =  1.  All  of  our  numerical  simulations 
assume  unit  magnetic  Prandtl  number,  t.e.,  p/(p0rj)  =  1,  so  we  subsequently  consider  only 
the  parameter  S  =  Sv  =  Sr. 

B.  Initial  Conditions 

The  initial  magnetic  and  velocity  fields  are  identical  to  the  vortex  system  of  Orszag 
and  Tang  (1979): 


and 


/>(x,  y,t  =  0)  =  constant  =  p°,  (6a) 


v(x, y,t  =  0)  =  —  siny  e*  +  sinx  ey,  (66) 


B(x,  y,t  =  0)  =  -siny  e*  +  sin2x  ey,  (6c) 
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where  er  and  ey  are  unit  vectors  in  the  x  and  y  directions  respectively.  For  reasons  given 
below,  our  investigations  of  emerging  shocks  use  an  initial  mass  density  which  can  differ 
from  unity.  The  initial  cross  helicity  and  normalized  correlation  between  the  velocity  and 
magnetic  fields  then  equal  0.5.  As  in  our  previous  calculations,  the  initial  mechanical 
pressure,  p(t  =  0),  comprises  a  mean  part,  ( p(t  =  0))  =  p°,  and  a  fluctuating  part, 
6p(x,y,t  =  0)  =  6p°(x,y).  Here  the  bracket  notation  represents  a  volume  average,  t.e.,  an 
average  of  a  function  (/)  over  the  entire  grid  (Q): 

f  f(x,y,t)  dx  dy 

G _ 

J  dx  dy 
G 

The  equation  for  the  dimensionless  speed  of  sound, 


(8) 


constrains  the  total  pressure  to  equal  or  exceed  zero.  The  lower  bound  on  p°  is  then  a 
positive  definite  functional  of  Sp°(x,y),  given  by  p^jn  >  max{—  6p0},  assuming  that  the 
choice  of  Sp°  is  not  identically  zero  in  the  computational  domain. 


Equation  (8)  indicates  that  formula  for  the  local  Mach  number  M  is 


=  . 

7  P 

Notice  that  the  pressure  must  be  greater  them  zero  everywhere  in  order  for  M  to  be 
real  and  free  of  singularities  throughout  the  grid.  To  obtain  a  “characteristic”  or  an 
“average”  Mach  number  M  at  time  t  =  0  for  the  flows  in  our  calculations,  we  substitute 
p  —  p°,  p  =  p°,  and  |v°j  =  >/(v2(t  =  0))  into  Eq.(9a): 

=  ^  . 

7  P° 

The  ratio  of  thermal  to  magnetic  pressure  depends  on  local  Mach  number  according  to 
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(10) 


3  =  =  2^lvl2 

\B\2  7A^2|B|2  ' 

In  all  calculations,  the  initial  “characteristic”  value,  /?c,  is  10/3,  to  permit  comparison  with 
the  M  =  0.6  case  considered  previously  [1].  The  quantity  fic  relates  to  0  in  Eq.(10)  as  M 
does  to  M.. 

Equation  (9a)  shows  that  one  could  obtain  flows  with  supersonic  regions  merely  by 
setting  the  initial  background  pressure  p°  close  enough  in  magnitude  to  the  minimum 
(negative)  value  of  the  fluctuating  pressure  6p°.  As  in  our  subsonic  calculations,  6p° 
would  be  the  solution  of  the  incompressible  MHD  equations  with  p°  —  1.  Unfortunately  as 
p°  — »  max{— <$p0},  the  supersonic  regions  move  inward  toward  the  center  of  the  grid,  which 
is  a  stagnation  point  of  the  flow  and  an  X-point  of  the  magnetic  field.  Figure  1(a)  shows 
the  initial  local  Mach  number  M  when  the  minimum  value  of  the  thermal  pressure  is  10-4. 
The  average  background  pressure  is  p°  «  2.6,  making  M  «  0.679.  The  initial  maximum 
value  of  M  is  approximately  2.2,  since  the  minimum  pressure  occurs  near  the  stagnation 
point.  Calculation  of  this  case  up  to  t  =  8  shows  only  small  differences  from  the  M  —  0.6 
subsonic  case  studied  previously,  even  though  the  latter  had  an  initial  maximum  value  of 
0.97  for  Ad.  In  fact,  both  calculations  exhibit  supersonic  flows  early,  but  no  shocks  occur 
[1].  Erlebacher  et  al.  have  reported  a  similar  result  for  3-D,  compressible,  Navier-Stokes 
turbulence  [14]. 

For  the  study  of  effects  related  to  emerging  shocks  within  this  model,  another  prescrip¬ 
tion  is  necessary.  Our  approach  is  to  increase  p°  above  1.0  while  retaining  the  definitions  of 
v°  and  B°  in  Eqs.(6b)  and  (6c).  In  addition,  we  allow  the  fluctuating  pressure,  6p° ,  to  sat¬ 
isfy  the  incompressible  equation  of  motion  with  p°  set  equal  to  one.  With  the  proper  choice 
of  p°,  we  can  define  a  set  of  initial  conditions  which  are  the  same  as  those  of  the  previous 
subsonic  study  except  for  the  initial  mass  density.  By  increasing  p°,  we  can  raise  the  char¬ 
acteristic  Mach  number  (Eq.(9b))  of  the  flow  above  1.0.  Since  these  choices  of  p°  and  6p° 
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do  not  satisfy  the  incompressible  equation  of  motion,  an  imbalance  exists  between  the  mo¬ 
mentum  gradients  of  the  fluid  and  the  opposing  magnetic  and  thermal  pressure  fields.  This 
imbalance,  in  fact,  causes  the  supersonic  disturbances  which  we  desire  to  study.  In  support 
of  this  approach,  Appendix  A  shows  that  attempting  to  maintain  a  balance  among  local 
forces  confines  us  to  subsonic  characteristic  Mach  numbers  M  <  0.679,  if  the  amplitudes 
of  the  velocity  and  magnetic  field  fluctuations  are  equal,  as  in  Eq.(6). 

Our  choice  of  initial  average  pressure  p°  is  10/3,  and  the  corresponding  characteristic 
flc  is  also  10/3.  These  values  are  the  same  as  in  the  M  =  0.6  subsonic  calculation  [1],  in 
which  p°  =  1.  The  present  calculations,  therefore,  are  linked  directly  to  our  previous  work. 
To  ensure  that  significant  supersonic  regions  are  present,  we  set  the  characteristic  Mach 
number  M  >  1,  giving  an  initial  mass  density  p°  >  2.78  by  Eq.(9b).  Here  we  calculate 
the  system  evolution  for  M  =  1.0,  p°  =  2.78,  and  M  =  1.5,  p°  =  6.25.  Figure  1(b)  shows 
the  initial  form  of  the  local  Mach  number,  «M,  for  the  M  —  1.0  calculation.  The  local 
Mach  number  varies  widely  over  the  system  from  a  value  of  0.0  at  the  stagnation  point 
in  the  center  of  the  computational  grid  to  a  value  of  1.6  within  sizable  supersonic  regions. 
Section  III  shows  that  these  regions  produce  shocks  which  alter  the  evolution  significantly 
from  the  subsonic  and  incompressible  baselines. 

C.  Numerical  Algorithm 

Previous  papers  have  described  our  numerical  algorithm  in  detail  [1,  15].  Our  code, 
CRXJNCH2D ,  implements  a  Fourier  collocation  method  and  employs  an  isotropic  trunca¬ 
tion  in  Fourier  space  at  each  time- level  [14].  The  modified  Euler  method,  a  second-order 
Runge-Kutta  scheme  [16],  discretizes  time.  The  time  step,  At,  is  limited  by  a  compressible 
MHD  Courant-Friedrichs-Lewy  (CFL)  number, 

Hen  =  SUp[M±£l+^]  <  o.3  ,  (11) 

where  Cs  =  sound  speed  and  h  =  2n/N  =  Ax  =  Ay.  When  shocks  are  present,  reduction 
in  the  time  step  according  to  Eq.(ll)  (as  compared  to  subsonic  flows)  often  is  insufficient 
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to  guarantee  a  smooth  solution  in  the  vicinity  of  the  shock.  In  such  situations,  the  thermal 
pressure  and  the  mass  density  can  have  local  negative  values,  and  increased  viscous  and  re¬ 
sistive  dissipation  often  are  necessary.  When  the  diffusive  terms  are  considered  separately, 
numerical  instability  can  occur  due  to  enhanced  transport  in  the  vicinity  of  the  shock 
discontinuity.  We  have,  therefore,  added  stability  criteria  based  on  viscous  dissipation  [17] 
and  resistive  dissipation,  respectively: 


K  = 


K 


N2At 

^Svpm  in 


<  1 


N2At 
4  Sr 


(12) 


Our  simulations  consume  approximately  17  ps  per  time  step  per  grid  point  on  a  single 
processor  of  the  Numerical  Aerodynamic  Simulation  Program  (NAS)  Cray  Y-MP  computer 
at  the  NASA  Ames  Research  Center.  For  M  =  1  and  Lundquist  numbers  S  =  50,  we  use 
2562  collocation  points  while  for  M  =  1,  S  >  100,  and  M  =  1.5,  S  =  50,  we  use  5122 
collocation  points.  The  calculations  with  5122  collocation  points  require  approximately  7 
million  words  of  core  memory.  With  these  dissipation  levels  and  resolutions,  CRUNCH2D 
can  simulate  processes  with  moderate  strength  shocks  (Ad  <  3.0),  as  verified  by  a  number  of 
standard  tests,  e.g.,  the  rupturing  diaphragm  problem  [15]  and  the  numerical  simulations  in 
I.  Our  model  shocks  are  approximately  5  computational  cells  thick  when  S  =  200  and  N  = 
2562  [15]. 


9 


III.  RESULTS  AND  DISCUSSION 


As  in  our  previous  paper  [1],  our  investigation  of  compressibility  in  the  Orszag-Tang 
vortex  model  relies  primarily  on  such  diagnostics  as  time-dependent  area  integrals,  contour 
diagrams,  generalized  streamline  plots,  and  wavenumber  spectra  of  various  autocorrelations 
and  cross  correlations.  The  reader  should  refer  to  that  paper  for  the  corresponding  plots 
when  M  ranges  from  0.2  to  0.6.  Supersonic  compression  and  the  presence  of  shocks  pro¬ 
foundly  influence  energy  dissipation,  the  time  at  which  the  onset  of  strong  reconnection 
occurs  at  the  magnetic  X-point,  the  longevity  of  reconnection  and  the  associated  current 
sheet,  and  the  evolution  and  extent  of  the  correlation  between  the  magnetic  and  velocity 
fields.  The  wavenumber  spectra  and  the  evolution  of  global  averages  thus  differ  signifi¬ 
cantly  from  those  of  subsonic  flows.  To  aid  in  the  interpretation  of  those  differences,  the 
discussion  begins  with  structural  features. 

A.  Structural  Effects  of  Supersonic  Flows  and  Emerging  Shocks 

For  the  incompressible  and  subsonic  systems,  the  structures  of  interest  are  those  as¬ 
sociated  with  magnetic  field  reconnection:  magnetic  islands  and  their  coalescence,  jets 
emerging  from  the  reconnection  zone,  vortex  quadrupoles,  and  electric  current  sheets  (“re¬ 
connection  current  sheets”  or  RCS).  In  the  supersonic  cases  discussed  below,  shocks  and 
their  associated  current  sheets  (“shock  current  sheets”  or  SCS)  and  distributions  of  vortic- 
ity  and  dilatation  (V  •  v)  constitute  additional  structural  features  which  become  dominant 
as  Mach  number  increases.  Since  the  system  is  decaying  with  time,  the  most  interesting 
structural  effects  occur  for  times  t  <  4.0.  Most  of  the  results  presented  in  this  subsection, 
therefore,  correspond  to  that  temporal  range. 

An  additional  process  emerges  due  to  the  self-interaction  of  the  magnetofluid  via  shock 
wave  propagation  through  the  system.  Because  of  the  symmetry  and  periodicity  of  the 
Orszag-Tang  model,  the  shocks  (and  later,  finite  acoustic  waves)  can  interact  multiple 
times  with  each  other  and  with  the  spatially  varying  mass  density  field.  The  latter  inter¬ 
actions  generate  progressively  smaller-scale  vortex  structures  and  density  fluctuations,  and 
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can  cause  the  magnetofluid  to  appear  turbulent  at  later  times  [1,  18,  19].  Picone  et  al.  [20] 
have  provided  a  simple  numerical  demonstration  of  this  phenomenon.  The  newly  generated 
vorticity  affects  the  local  mass  density  on  time  scales  corresponding  to  the  induced  rota¬ 
tional  motion  rather  than  the  shock  transit.  Thus  the  readily  detectable  manifestations 
appear  gradually  and  with  temporally  increasing  prominence. 

To  illustrate  the  shock-emergence  process,  we  follow  the  time  development  of  the  lo¬ 
cal  Mach  number  M.  throughout  the  grid.  Figure  2  shows  the  evolution  of  M.  for  the 
calculation  in  which  the  average  initial  Mach  number  M  =  1.0  and  S  =  100.  The  four 
regions  of  supersonic  flow  visible  in  Fig.  1(b)  have  moved  clockwise  and  outward,  com¬ 
pressing  the  fluid  ahead  of  them  and  forming  shocks  around  t  =  1.5.  Because  of  the 
periodic  boundaries,  identical  flows  squeeze  the  fluid  from  opposite  directions,  producing 
broad  regions  of  compression.  The  local  Mach  number  peaks  at  around  2.0  and  3.0  for 
M  ~  1.0,  S  =  100,  and  M  —  1.5,  S  =  50,  respectively.  The  shock  boundaries  expand, 
eventually  intersecting  at  the  central  X-point  (t  =  2.5),  compressing  the  weak  reconnec¬ 
tion  current  structure  residing  there  (Fig.  3),  and  producing  complex  transmitted  shock 
structures  (<  >  2.5).  The  transmitted  shocks  propagate  toward  the  boundaries,  interacting 
with  the  “parent”  shocks  and  weakening  through  dissipation  and  expansion. 

The  emerging  shocks  strongly  influence  the  spatial  variation  of  the  magnetic  field. 
Figure  3  shows  the  evolution  of  the  electric  current  density  over  the  interval  t  =  2  to  t  = 
4  for  M  =  1.0,  S  =  100.  Current  sheets  (SCS)  are  coincident  with  all  shocks;  the 
sharpest  SCS  reside  at  the  two  strongest  shocks.  Unlike  the  subsonic  and  incompressible 
systems,  therefore,  significant  dissipation  of  magnetic  energy  occurs  at  regions  other  than 
the  vicinities  of  magnetic  X-points.  The  magnetic  energy  available  to  drive  reconnection  at 
the  X-points  is  thus  lower  than  for  M  <  0.6.  At  t  =  2,  a  wide,  weak  reconnection  current 
sheet  (RCS)  is  present  in  the  center  of  the  grid,  where  strong  reconnection  and  sharp  RCS 
were  observed  much  earlier  in  the  subsonic  and  incompressible  simulations.  Compression 
of  this  central  region  occurs  when  the  shocks  intersect  at  t  =  2.5.  A  sharp  RCS  appears  at 
that  time,  indicating  that  significant  reconnection  and  resistive  heating  are  taking  place. 
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The  current  density  peaks  at  t  «  3,  later  than  in  the  subsonic  cases,  and  the  central  RCS 
is  prominent  through  t  =  6.  (Note  that  weak  current  sheets  also  appear  in  the  vicinity 
of  the  shocks  transmitted  through  the  center  of  the  grid.)  In  the  M  =  1.5  calculation, 
reconnection  is  delayed  more  significantly,  reaching  a  peak;  at  approximately  t  =  4. 

The  dilatation  evolves  similarly  to  the  shock  current  sheets  while  significant  vortex 
structures  reside  with  both  the  reconnection  current  sheet  and  the  shocks.  Figure  4  shows 
the  time  development  of  the  dilatation.  High  negative  values  near  the  shocks  indicate  that 
strong  compression  is  occurring  there.  The  weakening  compression  with  time,  evident  in 
the  decreasing  minimum  values,  is  due  to  the  relatively  low  Lundquist  numbers  of  the 
simulations  and  to  the  expansion  of  the  shocks  that  were  transmitted  through  the  central 
reconnection  region.  Vortex  structures  occur  in  regions  of  greatest  shock  curvature,  and 
an  additional  quadrupolar  vortex  resides  at  the  central  RCS,  as  in  the  incompressible 
and  subsonic  calculations  [1].  The  quadrupolar  vortex  structures  experience  increasing 
distortion  with  increasing  Mach  number.  Figure  5  shows  that,  when  the  shocks  reach  the 
center  of  the  grid  at  t  =  2.5,  shock  vortex  structures  diagonally  bisect  the  reconnection 
vortex  quadrupoles. 

In  the  incompressible  and  subsonic  systems,  the  alignment  (or  correlation)  of  the  ve¬ 
locity  and  magnetic  fields  grows  over  time  from  the  initial  value  of  50%  (cf.  I).  The  onset 
of  this  growth  occurs  slightly  later  in  the  subsonic  simulations,  the  delay  increasing  with 
M.  The  next  subsection  shows  that  the  supersonic  cases  undergo  similar  growth  of  cor¬ 
relation  after  much  larger  initial  delays.  The  period  of  growth  ends  with  a  plateau  or  a 
period  of  much  slower  growth.  The  correlation  value  at  the  plateau  is  lower  as  Mach  num¬ 
ber  increases.  Retardation  of  reconnection,  as  demonstrated  by  the  above  contour  plots, 
appears  to  be  a  significant  factor  in  reducing  the  rate  of  growth  and  level  of  correlation 
between  the  velocity  and  magnetic  fields  relative  to  the  cases  for  which  M  <  0.6.  Figures 
6(a)  and  (b)  show  streamline  and  magnetic  field  line  plots  at  t  =  3  for  M  —  1.0,  5  =  50. 
The  magnetic  field  has  been  compressed  along  the  horizontal  direction  and  stretched  in 
the  vertical  direction  relative  to  the  subsonic  systems.  The  shocks  themselves  cause  a  kink 
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in  the  magnetic  field,  but  their  presence  is  much  less  apparent  than  in  the  electric  current 
density,  which  depends  on  field  gradients.  The  field  lines  are  bent  toward  the  shock  fronts; 
hence,  we  are  observing  “fast”  shocks.  The  velocity  field  at  t  =  3  shows  the  early  stages  of 
reconnection  similar  to  the  situation  at  t  =  1  for  M  <  0.6.  In  particular,  an  X-point  forms 
in  the  velocity  field  near  the  center  of  the  grid,  where  the  X-point  in  the  magnetic  field 
resided  at  time  t  =  0.  This  X-point  in  the  velocity  field  remains  in  the  center  of  the  grid 
until  t  =  6.  By  that  time,  dissipation  has  significantly  reduced  both  the  intensity  of  recon¬ 
nection  and  the  effects  related  to  shocks.  In  Figs.  6(c)  and  (d),  the  velocity  and  magnetic 
fields  at  t  =  7  exhibit  mainly  a  large-scale  similarity.  We  address  this  more  thoroughly  in 
the  next  subsection  through  a  discussion  of  the  correlation  between  the  momentum  and 
magnetic  fields. 

The  explanation  for  the  delay  and  suppression  of  reconnection  lies  partly  in  the  out¬ 
ward  movement  of  the  supersonic  regions  due  to  the  initial  imbalance  between  the  cen¬ 
trifugal  foice  of  the  magnetofluid  and  the  opposing  thermal  and  magnetic  pressures.  This 
rarefies  the  region  near  the  central  magnetic  X-point  and  opposes  the  flow  of  fluid  into 
the  reconnection  region.  The  time  at  which  shocks  reach  the  center  of  the  grid  marks 
the  arrival  of  appreciable  fluid  from  the  outer  regions.  Only  then  does  one  observe  rapid 
reconnection  and  energy  dissipation,  as  seen  at  earlier  times  for  M  <  0.6.  The  presence  of 
current  sheets  at  the  shocks  also  weakens  reconnection  by  rapidly  dissipating  magnetic  en¬ 
ergy  that  would  otherwise  be  dissipated  at  the  magnetic  X-point.  The  transmitted  shocks 
propagating  away  from  the  reconnection  region,  although  weak,  impede  the  flow  of  fluid 
into  that  region,  somewhat  reducing  the  intensity  of  reconnection. 

In  summary,  our  simulations  show  that  the  emergence  of  shocks  outside  of  a  magnetic 
X-point  can  impede  reconnection  at  the  X-point  through  Ohmic  dissipation  of  magnetic 
energy  and  through  rarefaction  of  the  reconnection  region  itself.  Shocks  reaching  the 
X-point  compress  the  reconnection  region  and  signal  the  arrival  of  higher  density  mag¬ 
netofluid,  accelerating  the  reconnection  process.  We  expect  the  emergence  of  compressed 
regions  bounded  by  shocks  to  affect  reconnection  similarly  in  a  realistic  turbulent  MHD 
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flow  with  local  supersonic  regions  or  pressure-momentum  imbalances.  In  a  more  general 
situation,  reconnection  might  be  enhanced  in  one  region  by  compression  while  being  de¬ 
layed  by  rarefaction  or  shock-related  Ohmic  dissipation  in  a  neighboring  area. 

Figure  7  shows  the  evolution  of  the  mass  density.  Note  that  the  system  exhibits  180° 
rotational  symmetry,  as  does  the  initial  magnetic  field.  Fluid  and  magnetic  field  squeezed 
by  the  four  initial  supersonic  regions  (Figs.  1  and  2(a))  have  been  pushed  rapidly  along  a 
direction  tangential  to  the  bounding  shocks  (those  which  emerged  at  t  =  1.5)  and  inward 
toward  the  middle  of  the  system.  At  t  =  2,  the  shocks  at  the  upper  left-hand  (ULH) 
and  lower  right-hand  (LRH)  regions  extend  to  the  central  region,  and  the  remaining  shock 
structures  appear  to  include  Mach  stems.  The  magnetofluid  along  that  diagonal  (ULH 
shock  -  central  region  -  LRH  shock)  is  rarefied  while  the  compression  regions  in  the  other 
quadrants  consist  of  two  narrow  plateaus  separated  by  a  narrow  region  of  somewhat  lower 
density.  At  t  =  2.5,  the  ULH  and  LRH  shocks  have  intersected  at  the  center,  compressing 
the  reconnection  current  sheet  and  the  magnetofluid.  The  circular  patterns  at  the  left 
and  right  ends  of  the  RCS  correspond  to  the  ends  of  the  fluid  jet.  In  Fig.  7(c)  ( t  =  3), 
the  corners  of  the  transmitted  shocks  and  intersections  of  the  primary  or  incident  shocks 
have  introduced  additional  local  density  peaks  in  the  region  of  the  jet.  By  t  =  4,  the 
density  extrema  have  coalesced  and  spread  and  the  transmitted  shocks  have  expanded  and 
weakened,  interacting  with  the  primary  shocks  that  preceded  reconnection.  Figure  7(e) 
shows  the  mass  density  at  the  end  of  the  calculation  (t  =  8).  The  maximum  and  minimum 
values  indicate  that  the  mass  density  is  evolving  toward  a  uniform  state,  consistent  with 
the  dissipative  nature  of  the  medium.  However,  because  of  multiple  shock  interactions  with 
the  density  field  during  the  system  evolution  [20]  and  because  of  the  injection  of  small  scale 
fluctuations  by  magnetic  reconnection,  the  fluid  has  a  turbulent  appearance  with  sizable 
fluctuations  at  a  range  of  scales,  as  mentioned  at  the  beginning  of  this  section. 
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B.  Evolution  of  Global  Properties 
B.l.  Mach  Number  Dependence 

The  M  =  0.6,  5  =  50,  calculation  described  by  Dahlburg  and  Picone  [1]  provides  a 
useful  baseline  for  evaluating  the  evolution  of  global  properties  related  to  dynamic  align¬ 
ment,  to  production  of  vorticity  and  electric  currents,  and  to  energy  transfer  among  MHD 
variables  due  to  reconnection,  dynamo  action,  and  shocks.  Figure  8  shows  the  evolution 
of  severed  global  parameters  for  the  simulations  with  5  =  50  and  M  —  0.6,  1.0,  and  1.5. 

Figure  8(a)  shows  the  compressible  coefficient  of  correlation  between  the  momentum 
and  magnetic  fields: 


(pv  '  B) 

v'KpM)2)  (IBP)  • 
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For  the  compressible,  supersonic,  Orszag-Tang  vortex  system  studied  here,  ac  =  0.5  at 
t  =  0.  As  indicated  previously  [1],  the  growth  of  ac  relates  to  the  process  of  “dynamic 
alignment”  in  turbulent,  compressible  magnetofluids.  In  incompressible  models,  dynamic 
alignment  occurs  as  the  kinetic  and  magnetic  energy  decay  selectively  with  respect  to  the 
cross  helicity.  The  measure  often  used  in  incompressible  calculations  is,  therefore,  the  ratio 
of  the  cross  helicity  ( Hc )  to  incompressible  total  energy  (Ei),  given  by 


Hc  _  2(v  •  B) 

Ei  (MJ)  +  (IBP)  ' 
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In  all  of  our  investigations,  the  two  measures  of  correlation  have  qualitatively  similar 
evolutions  and  provide  useful  measures  of  differences  among  incompressible,  subsonic,  and 
supersonic  MHD  flows. 

Figure  8(a)  shows  that  the  baseline  subsonic  system  ( M  =  0.6)  undergoes  more  rapid 
growth  in  ac  than  do  the  supersonic  cases  (M  >1)).  In  fact,  the  M  =  0.6  calculation  has 
two  stages  of  rapid  growth  (for  t  <  8),  beginning  at  approximately  t  =  1.5  and  4.5,  respec¬ 
tively,  and  ending  with  plateaus  at  successively  higher  correlation  values.  The  correlation 
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q c  in  each  of  the  supersonic  calculations  exhibits  only  a  single  stage  of  rapid  growth  to  a 
plateau  for  t  <  8.  beginning  at  f  a  3.5  when  M  =  1.0  and  at  t  «  4.5  when  M  =  1.5.  The 
supe.  sonic  systems  also  experience  an  early  stage  in  which  the  correlation  level  wanders 
around  the  initial  value  of  0.5  before  the  onset  of  rapid  growth.  The  end  of  this  initial 
"induction  period"  occurs  at  roughly  the  same  time  as  the  peak  in  the  current  density  in 
the  reconnection  region  (Section  III. A.).  Incidentally,  the  M  =  0.0  curve  shows  an  initial 
period  of  slow  growth  which  also  ends  when  the  current  density  peaks  in  the  reconnection 
region,  in  this  instance  at  t  ss  2.0.  Figure  8(a)  shows  that  the  plateau  value  of  ac  at  late 
times  varies  approximately  linearly  with  the  reciprocal  of  the  characteristic  Mach  number 
M. 

In  the  compressible  Orszag-Tang  vortex  model,  growth  of  correlation  thus  appears 
to  depend  strongly  upon  the  dynamics  of.  and  energy  dissipation  by,  reconnection  at  the 
central  magnetic  X-point.  As  indicated  in  the  previous  subsection,  the  supersonic  flows 
in  the  Orszag-Tang  vortex  initially  rarefy  the  plasma  in  the  central  region,  weakening 
reconnection  there.  Shocks  bound  the  regions  of  high  density,  so  that  the  arrival  of  shocks 
in  the  reconnection  region  is  a  precursor  to  considerable  inflow  of  magnetofluid.  The  latter 
inflow  is,  in  turn,  necessary  for  significant  reconnection  to  occur,  and  ac  begins  sustained 
growth  shortly  thereafter.  The  increased  delay  in  sustained  growth  of  correlation  with 
increasing  Mach  number  thus  corresponds  to  the  later  arrival  of  the  shocks  at  the  central 
X-pcint.  In  addition,  as  the  characteristic  Mach  number  increases,  the  shocks  become 
stronger,  as  do  the  current  sheets  associated  with  them.  Then  the  amount  of  magnetic 
and  kinetic  energy  transformed  into  internal  energy  at  the  shocks  also  must  increase  with 
M .  This  leaves  less  energy  to  drive  reconnection  and,  consequently,  dynamic  alignment, 
causing  the  value  of  ac  at  a  given  time  to  be  lower  as  the  Mach  number  increases.  Following 
the  rapid  growth  of  correlation,  all  cases  reach  a  stage  of  very  slow  growth,  after  the 
reconnection  event  has  run  its  course. 

Figures  8(b)  and  (c)  show  the  temporal  evolution  of  the  kinetic  energy  (normalized  to 
one)  and  magnetic  energy  for  various  Mach  numbers.  The  dependence  of  both  on  M  >  1 
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is  much  stronger  than  for  the  subsonic  flows  (M  <  0.6),  which  showed  little  variation  with 
Mach  number.  In  incompressible  numerical  simulations  [10],  the  oscillations  in  kinetic  and 
magnetic  energy  were  approximately  180°  out  of  phase  and  were  interpreted  as  Alfvenic. 
Our  subsonic  calculations  [1]  were  quite  similar  to  the  incompressible  calculations  in  this 
regard,  and  the  same  interpretation  applied,  since  the  initial  kinetic  and  magnetic  energies 
were  equal  and  smaller  than  the  internal  energy  by  a  factor  of  at  least  five  and  because 
the  effects  of  acoustic  waves  were  small.  In  the  present  calculations,  however,  the  initial 
kinetic  energy  increases  as  M2  while  the  initial  magnetic  energy  remains  at  unity.  In 
addition,  shocks  and  magnetoacoustic  waves  are  present,  so  that  the  temporal  evolutions 
of  the  magnetic  and  kinetic  energies  are  not  linked  in  the  same  manner  as  for  the  subsonic 
and  incompressible  systems.  The  increased  magnitude  of  the  kinetic  energy  as  M  increases 
explains  the  greater  relative  sensitivity  of  the  magnetic  energy  to  supersonic  Mach  number 
(Fig.  S(c)  vs.  8(b)).  Small  percentage  changes  in  kinetic  energy  could  result  in  larger 
relative  changes  in  the  magnetic  energy,  whereas  small  fluctuations  in  the  magnetic  energy 
would  affect  the  kinetic  energy  even  less. 

Additional  phenomena  should  affect  the  evolution  of  the  magnetic  energy  shown  in 
Fig.  8(c).  First,  appreciable  reconnection  is  delayed  and  apparently  weakened  as  the  Mach 
number  increases.  Second,  given  that  we  axe  observing  fast  shocks,  compression  should 
enhance  the  magnetic  energy,  leading  to  larger  peak  values  •with  increasing  Mach  number. 
Because  of  the  larger  jump  in  tangential  magnetic  field  across  the  shock  with  higher  Mach 
number,  the  dissipation  of  magnetic  energy  by  the  shocks  and  their  associated  shock  current 
sheets  also  should  increase  with  Mach  number.  In  Fig.  8(c),  the  peak  magnetic  energy 
and  the  subsequent  dip  are  both  proportionately  larger  for  M  =  1.5  than  for  M  =  1.0. 
In  addition,  the  final  stage  of  decreasing  magnetic  energy  occurs  later  because  the  stage 
of  strong  reconnection  is  delayed.  Notice  that  a  shoulder  appears  at  t  «  3.5  for  M  =  1.0 
and  at  f  ~  4.5  for  M  =  1.5.  These  are  the  respective  times  at  which  the  rapid  growth  of 
correlation  began. 
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In  I,  our  interpretation  of  the  temporal  evolution  of  the  M  =  0.6  calculation  was 
different  from  that  in  the  previous  paragraph  because  no  shocks  were  present.  In  Fig. 
8(c),  the  magnetic  energy  is  enhanced  at  t  fa  2  when  M  =  0.6.  This  is  later  than  for 
the  supersonic  cases  and  is  consistent  with  an  Alfvenic  exchange  of  energy  from  the  mo¬ 
mentum  field  to  the  magnetic  field,  similar  to  that  occurring  in  incompressible  media  [1, 
21].  The  subsequent  decay  of  magnetic  energy  is  more  immediate  than  in  the  supersonic 
calculations,  and  no  “shoulder”  appears  in  the  curve.  We  attribute  this  to  the  fact  that 
strong  reconnection  is  not  delayed  when  M  is  subsonic  [1].  An  Alfvenic  transfer  of  energy 
between  the  momentum  and  magnetic  fields  is  likely  to  occur  in  the  supersonic  systems  as 
well,  and  we  cannot  state  with  certainty  that  compressional  effects  are  more  important. 
However,  the  stronger  Mach  number  dependence  of  the  supersonic  cases  leads  us  to  believe 
that  the  interpretations  in  the  previous  paragraph  are  reasonable. 

Figures  8(d)  shows  the  evolution  of  the  kinetic  enstrophy,  (u;2),  where  u  is  the  vor- 
ticity.  The  decrease  and  delay  in  kinetic  enstrophy  production  relative  to  the  subsonic 
calculation  (M  =  0.6)  must  relate  to  shock  formation  and  the  delay  in  appreciable  recon¬ 
nection,  as  discussed  previously.  Here  the  initial  rise  in  enstrophy  occurs  when  shocks  first 
emerge,  indicating  that  shock  curvature  is  contributing  significantly  to  vorticity  produc¬ 
tion,  according  to  Crocco’s  Theorem  [22].  The  structure  in  the  supersonic  systems  depends 
in  a  complicated  manner  upon  shock  formation  and  interactions,  as  well  as  upon  vorticity 
production  at  the  central  reconnection  current  sheet  and  the  secondary  RCS  residing  on 
the  boundaries. 

Figure  8(e)  shows  the  magnetic  enstrophy,  (j2),  where  j  is  the  electric  current  density. 
The  effects  of  shock  emergence  on  magnetic  energy  dissipation  appear  prominently  in 
the  form  of  a  massive  enhancement  at  t  «  1.5  for  M  =  1.5  and  as  a  smaller  peak  for 
M  =  1.0.  The  supersonic  cases  have  similarly  sized  peaks  at  t  «  3  ( M  =  1.0)  and 
t  «  4  ( M  =  1.5),  corresponding  to  reconnection  at  the  central  magnetic  X-point.  In 
contrast,  the  broad  enhancement  in  the  subsonic  system  indicates  that  steady  reconnection 
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occurs  at  approximately  the  same  level  over  the  interval  1.0  <t  <  3.0.  The  final  stage  of 
decay  occurs  later  as  the  Mach  number  increases. 

B.2.  Lundquist  Number  Dependence 

Figure  9  shows  global  averages  corresponding  to  M  =  1.0  and  5  =  50,  100,  and  200 
for  0  <  t  <  8.0.  In  Fig.  9(a),  the  compressible  correlation  coefficient,  ac,  shows  the  same 
qualitative  behavior  as  in  the  analogous  subsonic  and  incompressible  simulations  [1],  with 
the  degree  of  alignment  at  a  given  time  decreasing  as  S  increases.  The  kinetic  energy  (Fig. 
9(b))  also  displays  similar  behavior  to  its  subsonic  counterparts,  with  almost  no  Lundquist 
number  dependence  during  the  early  stage  of  development.  However,  the  timing  of  the 
fluctuations  in  average  kinetic  energy  is  different  from  that  of  the  subsonic  calculations 
(Fig.  5(c),  I).  The  fact  that  the  initial  kinetic  energy  in  Fig.  9  is  nearly  three  times  the 
initial  magnetic  energy  means  that  the  latter  affects  the  kinetic  energy  much  less  than  in 
the  subsonic  flows.  Apparently  the  compression  of  the  fluid  downstream  of  the  initial  high 
Mach  number  regions  and  the  formation  of  shocks  reduce  the  kinetic  energy  and  transfer 
energy  to  the  magnetic  field.  This  energy  transfer  occurs  through  Alfven  waves  in  an 
incompressible  magnetofluid. 

The  evolution  of  magnetic  energy  as  a  function  of  Lundquist  number  (Fig.  9(c)) 
is  less  similar  to  that  in  the  subsonic  cases.  The  initial  extremum  in  magnetic  energy 
occurs  earlier,  at  t  ss  1.5,  when  shocks  form,  as  opposed  to  t  «  2  for  M  <  0.6.  For 
S  =  100  and  200,  a  second  extremum  appears  at  t  «  4  which  is  not  seen  for  S  =  50  or  for 
M  <  0.6.  This  is  most  likely  due  to  a  combination  of  two  factors: 

(1)  Following  the  arrival  of  the  “incident”  shocks  at  the  central  magnetic  X-point  and  the 
onset  of  strong  reconnection  there,  transmitted  and  diffracted  fast  shock  waves  interact 
with  the  reconnection  current  sheet  and  with  remnants  of  the  incident  waves.  The  array 
of  shocks  no  longer  forms  a  network  bounding  quasistationary  regions,  as  seen  at  earlier 
times.  Instead  the  shocks  expand  and  propagate  through  a  large  portion  of  the  system, 
compressing  the  magnetic  field  and  increasing  the  magnetic  energy.  Interactions  among  the 
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various  shocks  and  magnetoacoustic  waves  enhance  the  complexity  of  the  wave  structure 
over  time.  The  coverage  of  the  grid  by  these  waves  thus  increases  with  time  for  t  >  1.5,  and 
their  influence  later  diminishes  as  dissipative  effects  become  more  prominent  for  t  >  4.0. 
This  sequence  is  implicit  in  the  broad  peak  observed  in  Fig.  9(b). 

(2)  The  reduced  rate  of  dissipation  of  shock  waves  and  other  magnetoacoustic  waves  at 
higher  Lundquist  numbers  increases  the  strength  of  the  compression.  The  broad  peak  in 
magnetic  energy  at  t  «  4  thus  becomes  more  prominent  as  Lundquist  number  increases. 

Figures  9(d)  and  (e)  show  that  the  kinetic  and  magnetic  enstrophies  increase  approxi¬ 
mately  linearly  with  Lundquist  number  and  have  similar  structures.  Distinct  peaks  emerge 
at  t  «  3  and  t  «  4  as  dissipation  decreases.  The  peak  at  t  =  3  corresponds  to  the  time  at 
which  the  current  in  the  reconnection  region  reaches  a  maximum.  As  in  Fig.  9(c),  the  peak 
in  current  density  at  t  =  4  reflects  the  wide  distribution  and  interactions  of  numerous  mag¬ 
netoacoustic  waves  and  shock  waves,  following  the  arrival  of  the  “incident”  shock  waves 
at  the  central  magnetic  X-point.  The  subsequent  onset  of  strong  reconnection  and  the 
transmission  and  diffraction  of  the  incident  shocks  through  the  reconnection  region  (e.g., 
Fig.  2(d))  results  in  a  complex  global  magnetic-field  structure  with  significant  gradients 
throughout  the  system  at  t  «  4.  Figure  3(d)  verifies  this,  showing  the  presence  of  sizable 
currents  over  a  large  portion  of  the  grid  at  that  time. 

C.  Spectral  Decompositions  of  Correlations 

Here  we  examine  the  calculation  with  M  —  1.0,  S  =  100,  although  the  trends  iden¬ 
tified  are  consistent  across  all  of  the  simulations  under  discussion.  Prior  to  the  formation 
of  shocks  at  t  ss  1.5,  the  wavenumber  spectra  are  typical  of  decay  problems  and  resemble 
those  of  the  subsonic  cases  [1].  Figure  10(a),  for  example,  shows  the  mass  density  auto¬ 
correlation  spectrum  at  t  —  1.  However,  after  shocks  form,  the  wavenumber  spectra  of  the 
mass  density,  velocity,  nonsolenoidal  and  solenoidal  components  of  the  velocity,  and  the 
magnetic  field  appear  quite  changed,  as  shown  in  Fig.  10(b)  and  (d)  -  (g)  for  time  t  =  2.0. 
Except  at  wavenumber  extremes,  each  is  approximately  linear,  with  a  slope  similar  to  that 
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of  the  reference  line  (A:-4).  The  magnetic  field  deviates  from  the  others,  showing  a  slightly 
steeper  slope.  The  magnetic  field  spectrum  also  differed  from  that  of  the  other  functions 
in  the  subsonic  simulations. 

This  change  of  spectrum  coincides  with  the  emergence  of  compressed  regions  bounded 
by  shocks.  With  the  passage  of  time,  the  high  wavenumber  portion  of  the  spectrum  falls 
from  the  reference  line  as  a  result  of  dissipation.  The  density  spectrum  in  Fig.  10(c)  exem¬ 
plifies  that  effect.  In  an  effort  to  model  the  spectral  features  of  these  compressed  regions,  J. 
Dahlburg  [23]  has  computed  numerically  the  equivalent  spectra  of  elevated  plateaus  with 
either  square  or  circular  cross  sections.  The  circular  enhancement  produced  an  approx¬ 
imate  k~ 3  autocorrelation  spectrum.  The  square  plateau  had  two  distinct  components: 
k~2,  characteristic  of  a  one-dimensional  planar  shock,  and  approximately  k~3  5.  Since 
the  actual  structures  have  cross  sections  similar  in  a  very  coarse  view  to  a  quadrilateral, 
the  square  “model”  is  closer  to  the  simulations.  These  considerations  show  that  simple 
idealizations  of  the  two-dimensional  compressed  regions  are  not  entirely  adequate  to  ex¬ 
plain  our  observations  quantitatively.  However,  the  qualitative  results  do  seem  to  explain 
the  change  in  autocorrelation  spectra  which  we  observe  when  two-dimensional  compressed 
regions  with  bounding  shocks  appear. 

Figure  11(a)  and  (b)  show  the  spectral  amplitude  of  the  cross-helicity  normalized  by 
the  spectral  amplitude  of  the  incompressible  energy  at  t  =  2.0  and  6.0,  respectively.  In  the 
subsonic  cases,  the  low  wavenumber  region  had  positive  values  while  the  high  wavenumber 
region  showed  an  anticorrelation,  although  this  bifurcation  was  less  pronounced  at  M  =  0.6 
than  at  M  =  0.2.  Here  we  see  that  the  velocity  and  magnetic  field  are  approximately  un¬ 
correlated  over  most  of  the  spectrum  at  times  prior  to  the  interaction  of  the  shocks  near 
the  X-point.  For  later  times  (t  >  4)  with  M  =  1.0,  the  lower  and  middle  wavenumbers 
show  a  positive  correlation  while  the  higher  wavenumber  components  are  approximately 
uncorrelated.  For  M  =  1.5,  5  =  50,  the  middle  and  higher  wavenumber  portions  of  the 
spectrum  are  uncorrelated,  while  the  lower  wavenumbers  have  positive  values  through¬ 
out  the  calculation.  The  presence  of  two-dimensional  compressed  structures  bounded  by 
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shocks,  the  coincident  dissipation  of  energy  by  the  shock  current  sheets,  and  the  delay  and 
reduction  of  energy  dissipation  by  reconnection  are  the  likely  causes  of  the  reduction  of 
correlation  at  the  smaller  scales.  In  addition,  these  factors  probably  are  responsible  for 
the  lower  globed  (average)  correlation  between  the  magnetic  and  velocity  fields  relative  to 
subsonic  flows. 

Mechanisms  which  suppress  the  growth  of  correlation  between  the  velocity  and  mag¬ 
netic  fields  at  temporal  scales  of  a  few  hours  or  less  axe  of  interest  in  interpreting  mea¬ 
surements  of  the  solar  wind  [24].  The  above  compressive  mechanism  can  occur  at  the 
scale  of  any  nonuniformity  within,  or  encountered  by,  the  solar  wind.  This  is  because  such 
nonuniformities  can  generate  shocks,  given  that  the  solar  wind  is  supersonic  and  is  not  a 
steady  flow.  Presently  the  excitation  of  Kelvin-Helmholtz  instability  by  velocity  shear  is 
the  leading  candidate  for  the  cause  of  a  reduction  in  correlation  between  the  magnetic  and 
velocity  fields  with  increasing  distance  from  the  Sim  [24].  We  point  out  that  the  interpre¬ 
tation  of  the  data  is  still  uncertain.  For  example,  early  data  analysis  indicated  that  MHD 
fluctuations  in  the  solar  wind  are  mainly  Alfvenic  and  propagate  away  from  the  Sun  [25]. 
Subsequent  theoretical  work  inferred  from  this  observation  that  the  population  of  either 
the  aligned  or  the  antialigned  state  would  become  dominant  with  increasing  distance  from 
the  Sun  (e.g.,  [26]).  Recent  data  analysis  has  indicated  a  much  smaller  population  of  out¬ 
wardly  propagating  Alfvenic  fluctuations  near  the  earth.  This  has  led  to  the  conclusion 
that  neither  state  of  alignment  dominates  at  distances  around  1  AU  from  the  Sun. 

The  correlation  of  the  fluctuations  in  thermal  and  magnetic  pressures  is  mostly  neg¬ 
ative  soon  after  the  shocks  have  formed,  as  exemplified  by  Fig.  11(c).  At  later  times 
( t  >  3.0),  this  effect  becomes  less  pronounced  or  disappears  entirely,  depending  on  M  and 
S.  According  to  arguments  by  Barnes  [27],  a  negative  correlation  is  consistent  with  the 
presence  of  quasistationary  structures  in  the  flow  field.  Clearly  such  structures  exist  here 
until  the  bounding  shocks  reach  the  central  magnetic  X-point,  after  which  time  the  shocks 
weaken  and  the  original  compressed  regions  change  significantly.  This  could  explain  the 
above  observation  that  the  negative  correlation  disappears  at  later  times.  In  their  studies 
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of  weakly  compressible  turbulence,  Shebalin  and  Montgomery  [6]  have  found  that  gradi¬ 
ents  of  the  total  pressure  (thermal  plus  magnetic)  relax  away  at  the  higher  wavenumbers. 
Their  studies  differ  from  the  present  ones  in  that  we  are  concerned  with  supersonic  flows, 
for  which  the  fluctuations  in  thermal  pressure  and  mass  density  are  not  necessarily  pro¬ 
portional. 

In  the  present  calculations,  the  mass  density  and  magnetic  pressure  showed  no  long- 
lived  or  readily  discernable  trends  in  spectral  correlation,  unlike  the  thermal  pressure  and 
the  magnetic  field.  We  can  understand  this  through  the  following: 

(1)  In  MHD  flows,  the  mass  density  changes  in  response  to  total  pressure,  which  includes 
contributions  from  both  magnetic  and  thermal  pressures,  and 

(2)  In  compressible  flows,  the  mass  density  and  thermal  pressure  fluctuations  need  not  be 
highly  correlated  due  to  the  influence  of  entropy  fluctuations. 

The  mass  density  and  magnetic  pressure  did,  however,  show  a  net  positive  (spatially  av¬ 
eraged)  correlation  during  the  evolution  of  the  system  at  average  Mach  number  M  =  1.5. 
At  such  high  Mach  numbers,  compression  of  the  magnetic  field  and  the  fluid  at  the  shocks 
most  likely  dominates  the  average  correlation  (over  the  grid)  and  could  be  sufficient  to  give 
a  net  positive  correlation. 
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IV.  Summary 

The  Orszag-Tang  vortex  system  provides  a  spatially  periodic  model  of  a  magnetofluid 
undergoing  transition  to  turbulence.  The  model  demonstrates  the  emergence  of  small 
scale  structure  through  magnetic  reconnection  and  current  sheet  formation.  Over  time, 
dynamic  alignment  of  the  momentum  and  magnetic  fields  occurs.  The  Orszag-Tang  sys¬ 
tem  also  represents  an  easily  reproducible  set  of  initial  conditions  for  comparisons  among 
different  numerical  models  of  MHD.  The  present  calculations  and  those  described  in  I 
extend  the  set  of  numerical  results  to  compressible  MHD  flows,  including  those  which  are 
moderately  supersonic.  We  have  found  that  the  presence  of  embedded  supersonic  flows 
and  the  existence  of  local  force  imbalances  at  moderately  high  0  can  cause  fast  shocks  and 
broad  compressions.  The  evolution  and  structure  of  the  “supersonic”  Orszag-Tang  system, 
therefore,  differ  significantly  from  those  of  its  subsonic  and  incompressible  counterparts. 
The  same  phenomena  will  occur  in  any  turbulent  magnetofluid  with  embedded  supersonic 
regions  and  local  properties  similar  to  the  Orszag-Tang  system.  As  a  consequence,  we 
predict  that  many  structural,  global,  and  spectral  properties  of  such  a  medium  will  differ 
substantially  from  those  of  subsonic  and  incompressible  turbulent  MHD  flows.  Here  we 
summarize  those  properties. 

In  the  Orszag-Tang  vortex  system,  the  presence  of  moderately  supersonic  regions 
(M  «  2)  can  result  in  the  emergence  of  shocks  in  the  region  outside  of  the  magnetic  X- 
point.  Accompanying  the  shocks  axe  prominent  current  sheets,  which  dissipate  magnetic 
energy.  Such  current  sheets  were  not  present  in  our  subsonic  and  incompressible  calcula¬ 
tions.  In  addition,  vortex  structures  reside  in  regions  of  changing  shock  curvature.  The 
shocks  occur  when  the  forces  d;.’e  to  thermal  and  magnetic  pressure  do  not  balance  the 
distribution  of  momentum  in  the  initial  conditions.  The  excess  centrifugal  force  pulls  fluid 
away  from  the  central  X-point,  rarefying  the  region  and  delaying  significant  reconnection 
there.  Since  the  shock  current  sheets  dissipate  magnetic  energy  at  early  times,  the  re¬ 
connection  is  weaker  as  the  initial  supersonic  Mach  number  increases.  The  intensity  of 
the  reconnection  increases  significantly  when  the  shocks  (and  the  dense  fluid  which  they 
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bound)  reach  the  region  of  the  original  X-point.  The  correlation  of  the  momentum  and 
magnetic  field  begins  growing  rapidly  thereafter.  Eventually  the  growth  slows  considerably 
so  that  the  correlation  reaches  a  “plateau,”  which  decreases  with  increasing  Mach  number. 
The  present  model  does  not  indicate  with  certainty  whether  reconnection  always  is  en¬ 
hanced  and  accelerated  by  compressive  fluid  motions  and  inhibited  by  rarefying  motions. 
However,  both  effects  occur  in  the  present  calculations. 

The  compressed  regions  bounded  by  shocks  dominate  autocorrelation  spectra  at  early 
times,  including  that  of  the  solenoidal  velocity  component.  We  have  not  yet  explained  why 
the  latter  occurs,  other  than  to  infer  that  the  solenoidal  and  nonsolenoidal  components  are 
closely  coupled  and  to  note  that  solenoidality  of  the  flow  does  not  imply  incompressibility  of 
the  medium.  As  opposed  to  the  subsonic  and  incompressible  cases,  the  normalized  spectral 
correlation  between  the  velocity  and  magnetic  fields  is  near  zero  over  the  middle  and 
upper  portions  of  the  wavenumber  domain.  The  correlation  between  thermal  and  magnetic 
pressure  is  negative  over  a  wide  wavenumber  range  at  early  times  for  all  calculations  with 
M  >  1.0.  This  is  consistent  with  the  presence  of  broad,  long-lived  regions  of  compression. 
At  later  times,  these  regions  are  disrupted  and  the  spectral  anticorrelation  of  thermal  and 
magnetic  pressures  disappears. 

To  complete  our  investigation  of  the  compressible  Orszag-Tang  model,  we  are  currently 
investigating  the  effect  of  varying  /?  on  the  above  conclusions.  The  present  and  previous 
simulations  set  /?  >  10/3,  so  that  the  magnetic  field  has  not  played  as  dominant  a  role  as 
it  might  in  some  space  or  laboratory  plasmas.  Our  next  step  in  the  study  of  compressible 
MHD  turbulence  will  be  to  investigate  the  decay  of  magnetofluids  from  random  initial 
conditions  with  and  without  embedded  supersonic  flows. 
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APPENDIX  A 

Pressure  Imbalance  in  the  Supersonic  Orszag-Tang  Vortex  System 

The  present  model  of  supersonic  regions  within  the  Orszag-Tang  vortex  system  uses 
a  fluctuating  pressure  p i  =  6p°  which  is  consistent  with  Eqs.(6)  for  an  incompressible 
magnetofluid  of  initial  mass  density  p°  =  1.  The  corresponding  velocity  and  magnetic 
fields,  denoted  by  Vj  and  Bi,  respectively,  then  satisfy  the  elliptic  equation: 

V.(v,  Vv1)  =  i[— i  V2(p,  +  IB, |2)  +  V-(B,.VB,)]  .  (A.l) 

P  & 

Equation  (A.l)  represents  the  force  balance  between  the  motions  of  the  magnetofluid  and 
the  confining  thermal  and  magnetic  pressures.  The  flat  background  pressure,  p°,  and  the 
characteristic  Mach  number,  M,  must  be  consistent  with  p°  and  |v°|  =  y/(vfj,  through 
the  equation 


M2  = 


2p°lv°l2 

7  P° 


(96) 


For  the  subsonic  cases  with  p°  =  1,  Eqs.(A.l),  (6),  and  (9b)  represent  a  consistent  speci¬ 
fication  of  the  initial  conditions.  As  stated  in  Section  II.B,  the  constraint  on  the  average 
pressure  p°  limits  its  minimum  value  to  max{—  pi(x,  y)}  +  e,  where  e  >  0,  so  that  the  local 
Mach  number  will  be  defined  everywhere.  For  7  =  5/3,  Fig.  1  of  Dahlburg  and  Picone 
[1]  shows  that  p°  is  >  2.6.  Equation  (9b)  then  limits  the  characteristic  Mach  number 
to  M  <  0.679.  Even  when  p°  is  sufficiently  close  to  2.6  for  supersonic  regions  to  exist 
within  the  initial  flow  field,  we  do  not  observe  shocks,  and  the  quantitative  behavior  of 
such  systems  is  quite  close  to  that  of  the  M  =  0.6  case  presented  in  I. 

Equation  (9b)  shows  that  one  could  increase  the  characteristic  Mach  number,  Af,  by 
increasing  the  amplitude  of  the  initial  velocity  field,  increasing  the  initial  mass  density 
or  lowering  p°.  If  the  initial  velocity  in  Eq.(6)  were  increased  by  a  constant  factor,  we 
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would  also  increase  the  magnetic  field  amplitude  by  the  same  factor  in  order  to  maintain 
the  initial  normalized  cross  helicity  in  Eq.(14)  at  0.5.  Given  this  constraint,  consider  each 
alternative  separately: 

(1)  Scale  the  velocity  field  up  by  a  constant  factor.  Set  B'  =  aBi  and  v'  =  ^Vi  with 
a  >  1.  Then  by  Eq.(A.l),  p'  =  1,  and  p'  =  a2pi  =*>  p0'  >  2.6a2.  By  Eq.(9b),  we  have 
M  <  0.679. 

(2)  Scale  the  mass  density  up  by  a  constant  factor,  so  that  p'  =  a  with  a  >  1.  Then  by 
Eq.(A.l),  a  solution  is  v'  =  "Vi/y/a,  B'  =  Bj,  and  p'  =  p\  =>  p0'  >  2.6.  By  Eq.(9b),  we 
have  M  <  0.679.  Note  that  this  solution  also  fails  to  maintain  the  initial  normalized  cross 
helicity  in  Eq.(14)  at  0.5. 

(3)  Scale  the  initial  pressure  down  by  a  constant  factor,  so  that  p'  =  api  and  p0'  >  2.6a, 

where  a  <  1.  Equation  (A.l)  gives  the  solution  v'  =  >/avi,  B'  =  and  p'  =  1. 

Again  Eq.(9b)  shows  that  M  <  0.679. 

This  shows  that  one  cannot  increase  M  to  supersonic  values  while  maintaining  the 
above  constraints  on  v  and  B  and  satisfying  Eq.(A.l).  We  can  achieve  supersonic  values 
of  M  only  by  violating  the  balance  of  forces  and  pressures  implied  by  Eq.(A.l).  Hence, 
our  approach  is  to  increase  p°  above  2.78  and,  at  the  same  time,  to  allow  6p° ,  v,  and  B 
to  satisfy  Eq.(A.l)  with  p°  set  equal  to  one.  Since  the  fluctuating  pressure  will  not  balance 
the  other  local  forces,  finite  magnetoacoustic  waves  or  shock  waves  will  emerge  in  the  flow. 
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Fig.  1  —  Contours  of  constant  local  Mach  number  Jt  at  time  t  =  0.  (a)  Characteristic  Mach  number  M  =  0.679  with 
balanced  forces  throughout  grid.  Minimum  value  of  the  initial  pressure  is  1  x  10~4.  Minimum  and  maximum  values  of 
JC  are  0.0  and  2  185.  (b)  Characteristic  Mach  number  M  =  1.0  with  unbalanced  forces  in  the  grid.  Minimum  and 
maximum  values  of  Jt  are  0.0  and  1.620. 


Fig.  2  -  Evolution  of  the  local  Mach  number  for  M  =  1.0  and  S  =100  (minimum,  maximum  values):  (a)  r  =  1.5  (0.0. 
1.970),  (b)  i  =  2.5  (0.0,  1.659),  (c)  t  =  3.0  (0.0,  1.932),  (d)  t  =  4.0  (0.0,  1.283). 
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g.  3  —  Evolution  of  the  electric  current  density  for  M  =  1.0  and  S  =  100  (minimum,  maximum  values):  (a)  t 
13.1,  21.6),  (b)  t  =  2.5  (-17.6,  28.0),  (c)  t  =  3.0  (-15.6,  38.1),  (d)  t  =  4.0  (-17.6,  22.3). 
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Fig.  6  —  Plots  of  the  generalized  stream  function  for  M  =  1.0  and  S  =50  giving  direction  of  velocity  ((a)  and  (c))  and 
magnetic  field  ((b)  and  (d))  at  respective  times  r  =  3.0  and  7.0. 
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Fig.  8  —  Evolution  of  global  averages  for  S  =50  and 
M  =  0.6  (A),  M  =  1.0  (O),  and  M  =  1.5  (□).  (a) 
Compressible  correlation  coefficient,  ac,  (b)  Kinetic 
energy,  normalized  to  one  at  r  =  0.  (c)  Magnetic 

energy,  (d)  Kinetic  enstrophy,  { oi 2).  (e)  Magnetic 
enstrophy,  U2). 


Fig.  9  —  Evolution  of  global  averages  for  M  =  1.0  and 
S  =50  (A)  S  =  100  (O),  and  S  =  200  (□).  (a) 
Compressible  correlation  coefficient  ac,  (b)  Kinetic 
energy,  (c)  Magnetic  energy,  (d)  Kinetic  enstrophy. 


Fig.  10  —  (Continued)  Autocorrelation  spectra  (log  of  magnitude  squared)  vs.  log  of  squared  wavenumber  for  (a)  Mass 
density,  t  =  1.0,  (b)  Mass  density,  t  =  2.0,  (c)  Mass  density,  /  =  4.0,  (d)  Velocity,  t  =  2.0,  W  Nonsolenoidal  component 
of  the  velocity,  t  =  2.0,  (f)  Solenoidal  component  of  the  velocity,  /  =  2.0,  (g)  magnetic  field,  i  =  2.0. 


